Asymptotic expansion of the inverse of sums of matrices

Theorem : For a matrix $C$ that can be expressed as $C = A + \varepsilon B$, the inverse matrix $C^{-1}$ can be expressed as

\begin{equation} C^{-1} = A^{-1} - \varepsilon A ^{-1} B A^{-1} + \mathcal{O}(\varepsilon^2) \end{equation}

Proof : The matrix $C$ obeys the Woodbury matrix identity,

\begin{equation} (A + \varepsilon B)^{-1} = A^{-1} - A ^{-1} (I + \varepsilon B A^{-1})^{-1} A^{-1} \end{equation}

Via the Neumann series, we get that

\begin{eqnarray} (I + \varepsilon B A^{-1})^{-1} &=& \sum_{k = 0}^\infty (-\varepsilon B A^{-1})^k &=& I - \varepsilon B A^{-1} + \sum_{k = 0}^\infty \varepsilon^{k + 2} (- B A^{-1})^{k + 2} \end{eqnarray}

The last term is indeed $\mathcal{O}(\varepsilon^2)$.